Disorder to order: how halide mixing in MAPbI3−xBrx perovskites restricts MA dynamics

Mixed-halide lead perovskites are of particular interest for the design of tandem solar cells currently reaching record efficiencies. While halide phase segregation upon illumination of mixed perovskites is extensively studied, the effect of halide disorder on A cation dynamics is not well understood, despite its importance for charge carrier diffusion and lifetime. Here, we study the methylammonium (MA) reorientational dynamics in mixed halide MAPbI3−xBrx perovskites by a combined approach of experimental solid-state NMR spectroscopy and molecular dynamics (MD) simulations based on machine-learning force-fields (MLFF). 207Pb NMR spectra indicate the halides are randomly distributed over their lattice positions, whereas PXRD measurements show that all mixed MAPbI3−xBrx samples are cubic. The experimental 14N spectra and 1H double-quantum (DQ) NMR data reveal anisotropic MA reorientations depending on the halide composition and thus associated disorder in the inorganic sublattice. MD calculations allow us to correlate these experimental results to restrictions of MA dynamics due to preferred MA orientations in their local Pb8I12−nBrn “cages”. Based on the experimental and simulated results, we develop a phenomenological model that correlates the 1H dipolar coupling and thus the MA dynamics with the local composition and reproduces the experimental data over the whole composition range. We show that the dominant interaction between the MA cations and the Pb–X lattice that influences the cation dynamics is the local electrostatic potential being inhomogeneous in mixed halide systems. As such, we generate a fundamental understanding of the predominant interaction between the MA cations and the inorganic sublattice, as well as MA dynamics in asymmetric halide coordinations.


Exact amounts of precursors for mechanochemical synthesis
MAPbI 3 : 1.900 g of MAI (11.95 mmol) and 5.509 g of PbI 2 (11.95 mmol) were milled for 50 min to prepare the black powder. MAPbI 2 Br 1 : 5.150 g MAPbI 3 (8.35 mmol) and 2.000 g MAPbBr 3 (4.18 mmol) were milled for 50 min to prepare the black powder.  Table S1. Experimental lattice constants for MAPbI 3-x Br x (x= 0, 1, 1.5, 2, 3) extracted by refinements of the PXRD patterns depicted in Fig. 1a. Lattice constant a / Å MAPbI 3 6.27 Modelled as pseudo-cubic MAPbI 2 Br 1 6.17 MAPbI 1.5 Br 1. 5 6.11 MAPbI 1 Br 2 6.04 MAPbBr 3 5.93  As the refocusing time of 10 µs used for the 207 Pb spin echo experiments (Fig 2 and Fig S2) is significantly shorter than the 207 Pb spin spin (T 2 ) relaxation determined for MAPbI 3 of ~40 µs 4 or in MAPbBr 3 (T 2 *=70 ms) 5 , we do not expect significant changes in the 207 Pb NMR spectra due to T 2 relaxation. Thus, we evaluate the 207 Pb spectra as pseudo-quantitative and compare the experimental integrals of PbI 6-x Br x environments (x=0, 1, …, 6) to the ones expected from the statistics for a random distribution of halides (Tab. S2-S4). This is supported by the fact that the halide composition obtained from XRD data matches those obtained from the NMR spectra for each halide composition.   Figure S6. 1

First-principles C Q calculations
The time-averaged EFG is calculated from the MD trajectory for x = 1.5. The EFG is calculated from the self-consistent potential using the method of Petrilli et al. 10 as implemented in VASP. The EFG tensor of each 14 N is averaged over the MLFF-MD trajectory and, after this averaging, rotated onto its principal axes to obtain C Q and η, using the quadrupole moment from the Pyykkö compilation. 11 The EFG tensors are calculated at intervals of 1 ps, in the 4 × 4 × 4 supercell, requiring a total of 355 calculations. In order to keep the computational load manageable we use the PBE functional 12 instead of the SCAN that was used for the ML potential. We used the PBE PAW potentials "Pb", "I", "Br", "C", "N", "H" from the VASP database, so a Pb potential with only 4 unfrozen (valence) electrons. As we do not need to optimize structures and focus on the N electronic structure, this is a reasonable approximation.

Machine Learning Force Field Method 5.1 MLFF: MAPbI 3-x Br x supercells
The MD simulations have been performed on 4x4x4 supercells containing 64 formula units as shown in Figures S8 and S9. Only the halide ordered systems are shown. The random systems have been generated out of the pure iodine systems, by randomly replacing iodine ions with bromine ions unit the desired fraction x is reached. Figure S8. The pure halide x=0 (left) and x=3 (right) supercells. Figure S9. The halide ordered x=1.5 supercells: structure A (left) and B (right).

MLFF: Calculation of average dipolar coupling and it's relation to second moment
We follow the approach as presented by Goc et. al. to numerically calculate the van Vleck second moment for the dipolar interactions of the protons. Experimentally the average dipolar interaction is extracted from double-quantum built-up curves following the approach by

MLFF: Fitting D intra and D inter
We have fitted the function for each of the seven systems as function of  In an attempt to characterize the halide distributions over these cages with just a few numbers, we introduce an order parameter O=[N xy ,N xz ,N yz ], containing the number of iodine ions (N=0...4) in each of the three cartesian planes through the center of the cage as defined in Figure 4c. This order parameter has been added to each of the 64 environments shown in Figures S13-S17. We further calculate the standard deviation of its three components, thus expressing the inhomogeneity and asymmetry of the halide coordination of a MA molecule in just a single number. We observe that in general the distributions corresponding to configurations that have a large (like [124]) show strongly preferred orientations, whereas configurations with a zero (like [222]) show a close to flat distribution of the C-N axis orientation, i.e., there is hardly any preferential orientation of the MA cations in such environments. Figs. S16-S18 explore these findings in more detail. They assess the correlation between and the corrugation of the polar distribution .

Figure S11. Average variance of the order parameter observed in the infinitely large supercell computed numerically. The entropy of mixing of a binary gas is a good approximate of this function.
We have numerically calculated the function in the limit of very large supercell sizes and ( ) plotted the result in Figure S11. Here, we applied the following (Python) routine: For every x we generate 10,000 random environments [XYZ], calculate the variance of their order parameter and average them. As shown in the figure, the mixing entropy of a binary gas multiplied by 1.045 is a good approximate of this numerical result. Note that for the mixing entropy formula the fraction expressed in the range x=[0..1] is used.

MLFF: Polar distributions of MA molecules
The orientation of each of the 64 molecules in the supercell is extracted as function of MD simulation time. Figure S12 illustrates the reference frame in which the orientation is In general, we observe that the MA molecules in the ordered structures show more highly preferred orientations as compared to the randomly distributed halide structures. This means that the molecules in these systems will rotate less isotropically and will therefore result in a higher D intra contribution to the 1 H-1 H dipolar coupling as compared to the random structures.
Note that the 64 distributions are ranked based on the level of anisotropy of the local cage that the particular molecule is in. This level is calculated by the standard deviation of the order parameter. For example, Figure S15 starts

MLFF: Dependence of the corrugation of the polar distributions on the order parameter
We will now assume that a large part of the symmetry and corrugation of the polar distribution corresponding to a particular molecule can be related to its local environment of the 12 nearest surrounding halides. We propose the following two physical parameters to model the corrugation level of the polar distribution: We have tested the correlation between these two parameters and the variance of the polar distribution. This variance is a direct measure for the corrugation of the polar distribution. As a correlation test the Pearson correlation coefficient is calculated. The resulting coefficients are summarized in Table S7. A value of -1 or 1 would indicate perfect (anti-)correlation. A coefficient of zero indicates no correlation, as in the case of two random sets of numbers. What the coefficients in Table S7 show is that the local concentration does not correlate with the level of corrugation. For the three different supercells different signs of the Pearson coefficients are observed. However, the level of anisotropy does correlate with the level of corrugation, and shows positive Pearson coefficients that increase with the global Bromide fraction. Note that these correlation coefficients are small, but statistically significant. It illustrates that the composition of the local cage is important, but the level of anisotropy is not fully sufficient to predict the polar distribution of a molecule.  ( ) The least-mean-square fitted red linear curve in Figures S18-S20 clearly shows that there is a correlation between the corrugation, and therefore the level of non-isotropy, of the polar distribution of the molecule depending on the local environment that it is in. Therefore, it is a legitimate descriptor to model the level of non-isotropy of the polar distribution. As expected, we also see that this simple descriptor (the local order parameter O) is not perfect. There is a large spread in the level of corrugation of polar distributions assigned to environments that have the same value. Note that O does not describe how the I and Br ions are distributed in the cartesian planes, nor does not take into account the type of the nearest neighbour environments. Of course, all these effects could be included, but it will greatly increase the complexity of the order parameter.

MLFF: Reorientation dynamics of MA molecules
The autocorrelation functions of the molecular orientations p(t) as sketched in Figure S12 are calculated for all molecules in the random structures. This is done in the same manner as explained in Ref. 16 . The autocorrelation functions (blue lines) are shown in Figures S22-S24 and fitted with a double exponential (orange lines). The fit is very good for all molecules in the x=1.125 and 1.5 random systems, and sufficient in the x=1.875 system. Distributions of all fit parameters are shown in Figure S21. In all cases we observe a fast 'thermal' decorrelation process of between 0.5-3 ps and a slow 'flipping' decorrelation process of between 3-30 ps.